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Abstract 

This paper is devoted to the approximation of a nonstandard Darcy problem, which 
modelizes the ow in porous media, by spectral methods: the pressure is assigned 
on a part of the boundary. We propose two variational formulations, as well as three 
spectral discretizations. The second discretization improves the approximation of the 
divergence-free condition, but the error estimate on the pressure is not optimal, while 
the third one leads to optimal error estimate with a divergence-free discrete solution, 
which is important for some applications. Next, their numerical analysis is performed 
in detail and we present some numerical experiments which conrm the interest of the 
third discretization. 
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Introduction 

We consider the following Darcy problem: 


u + Vp = f 

in SI, 

(1.1) 

div u = 0 

in Q, 

(1.2) 

£ 

II 

3 

on r\ 

(1.3) 

P = V 

on r 2 , 

(1.4) 


where Q is the plane square ] — 1,1[ 2 and n = ( 711 , 712 ) is the exterior unit normal to the 
boundary T = dQ. The boundary T is divided into two parts: the horizontal portion T 1 = 
{(x,y )| — 1 < x < l,y = ±1} and the vertical portion T 2 = {(x,y)\x = dzl, —1 < y < 1}. 
As we can see, the boundary condition on T 2 are nonstandard, since we prescribe the value 
of the pressure on T 2 . On the contrary, we have a classical condition on the portion T 1 . 
These conditions are described in the following figure. 

Boundary condition 
on the normal velocity 

T 1 


Q Boundary condition 

T 2 T 2 on the pressure 


T 1 

Boundary condition 
on the normal velocity 

Figure 1.1 

The equations of the Darcy problem not only modelize the flow in porous media, but 
also appear in the projection techniques for the solution of Navier-Stokes equations (see 
[10] and [15]). The nonstandard boundary conditions, where the pressure is assigned on 
a part of the boundary, have a physical meaning: typically the portion T 1 corresponds to 
rigid walls, whereas the entry or exit of the fluids takes place through T 2 . The spectral dis¬ 
cretization with this type of boundary conditions was only studied within the framework of 
the Stokes problem (see [6], [4] and [5]), while Azaiez, Bernardi and Grundmann proposed 
in [2] the spectral discretization of the standard Darcy problem, where the normal velocity 
is assigned on the boundary. 

This paper is devoted to the spectral discretization of the nonstandard Darcy problem. 
First, we give two variational formulations. Each one leads us to well-posed problems. Se¬ 
cond, we study the regularity of the solution by using a mixed problem of Dirichlet- 
Neumann for the Laplace operator. Next, from the first variational formulation, we derive a 
first spectral discretization, which is simple, out, in order to improve the approximation of 
the divergence-free condition, we study a second spectral discretization, where the inf-sup 
condition is obtained with more difficulty and where the error estimate on the pressure is 
not optimal. Finally, the second variational formulation yields a third spectral discretiza¬ 
tion, which leads us to optimal error estimate and a divergence-free discrete solution. 


Boundary condition 
on the pressure 
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An outline of this paper is as follows. The two continuous variational problems, as 
well as the regularity of the solution, are studied in Section 2. Section 3 is devoted to 
the analysis of three spectral approximations of this problem in the case of homogeneous 
boundary conditions. In Section 4, we present the algorithms that are used to solve the 
first and third discretizations, together with some numerical experiments. 


Statement of the problem and notation 

In order to set this problem into adequate spaces, recall the definition of the following 
standard Sobolev spaces (cf. J. Necas [13]). For any multi-index k = (fcijfo) with ki > 0, 
set \k\ = ki + k 2 and denote 

8 V = - t - r . 

dx^dx^ 

Then for any integer m > 0 and any plane domain £1 whose boundary is Lipschitz- 
continuous(cf. Grisvard [12]), we define: 

= {v € L 2 (fi); d k v € L 2 (Q) for 1 < \k\ < m} , 
equipped with the seminorm 

= ( £ L\d k v\ 2 (ix)5, 

\k\=m Jn 

and norm(for which it is an Hilbert space) 

m 

IM|if"(sJ) = (EE \\d kv \\h(n)) 1/2 - 

|fc|=0 k 

For extensions of this definition to non-integral values of m (see [11,12]), let s a real 
number such that s = m + cr with m £ IN and 0 < a < 1. We denote by H S (Q) the space 
of all distributions u defined in Q such that u E H m (Q) and, V|a| = m, 

f f (d a u(x) - d°u(y)) 2 . , 

LL \\x-y\\^ dxdy<+ °°- 

It can be shown that H s (£l) is a Hilbert space for the scalar product 

. . . , r r (d a u(x) — d a u(y))(d a v(x) — d Q v(y)) 

(U,v)„ a +E/ n / n l \\x y\\^ dxdy - (21) 

Let T' be an open part of the boundary dfl of class £7 m-1,1 and T[ the mapping 
v ^ V|r' defined on H m (Q). We denote by H m ~i( T') (see [7,12]) the space f2)) 

which is equipped with the norm: 

M H m-i (F') = inf{||i;|| ff m (n) , v G H m (Q) and vp* = ip}. (2.2) 

In this text, we shall use the spaces i7 1 / 2 (r / ) and H 3 / 2 (T') corresponding to m = 1 and 2. 

Let us define the space Hqq 2 (F') = {i?|r', v E i7 1 (Q), Vx E dQ \ T', ^|an(x) = 0}. We 
shall also be interested in the dual space of Hqq 2 (T'), 

H - 1/2 (r') = (-ffoo 2 (r'))'- (2.3) 

We shall use the Hilbert space i7(div;f2) = {v E L 2 (Q) 2 : div v E L 2 (fl)}, equipped 
with the norm 
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IM|if(div;Sl) — ((IM|l 2 (S 2)) 2 + (lldiv v|| i 2 ( S! )) 2 ) 2 . (2.4) 

For vanishing boundary values, we define: 

Hi(Q) = {veH 1 {Q) : V\ dn = 0}. 

For A =] — 1,1[, we denote the norm in L 2 (A) by |HIo,a — {J ( v ( x )) 2 dx)*, the semi- 
norm in if 1 (A) by Hi ,a = {J (v'(x)) 2 dx )2 and the norm in H 1 ( fi) by 
||^ ||i, a = (J i ((v(a:)) 2 + (v'(x)) 2 )dx)*. 

We note x = (x, y) the generic point of the square Q and we call T/, T//, Tjjj 
and Tjv the edges of Q, starting from west and turning counterclockwise. For each J, 
J = I, II, III, IV, the extremities of the edge Tj are aj_i and aj, with the convention 
a 0 = aiv, the exterior unit normal vector to Tj is denoted by nj and the counterclockwise 
unit tangent vector is tj. Figure 1.2 below presents this notation. 

For any domain A in IR or IR 2 and for any nonnegative integer n, P n (A) stands for 
the space of all polynomials on A with degree < n with respect to each variable. We also 
use the notation P°(A) for the subspace P n (A) n Hq(A). For A =] — 1,1[, the family 
(L n ) n of Legendre polynomials is a basis of the spaces P(A) of polynomials on A (we refer 
to [7, Chap. I] for the properties of the orthogonal polynomials). These polynomials are 
orthogonal to each other in L 2 ( A) and are caracterized as follows: for any integer n > 0, 
the polynomial L n is of degree n and satisfies L n (l) = 1. Let us recall some properties 
that we need. The family ( L n ) n is given by the recursion relation: 


Lo = V Li(C) = C, 

(n + 1)L„ +1 (C) = (2n + 1)C-MC) - nL n _ i{C), n> 1. 


Each polynomial is a solution of the differential equation 

((1 — ( 2 )L' n )' + n(n + l)L n = 0, n > 0, 
and its norm is given by 

l|inll “' A “2nTl’ n - 0 ' 


ao — aiv 


ni 


niv 


ai 


\',v T IV 

r ui 

Tl 


0 

Till 

T/ 


TU 

r n 


nii 


am 


nm 


an 


Figure 1.2 

Three consecutive polynomials are linked by the integral equation 


(2.5) 


( 2 . 6 ) 

(2.7) 
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f( L„(f) d£ = - L„_i(0), 

From (2.6) and integration by parts, we derive 

L' n ( 1) = " (! ^, \\L' n \\l ,A = n(n + l), 

6 P„(A), \<Pn\lA < ^ / 3 « 2 ||<Pn||o,A- 


n > 1. 


( 2 . 8 ) 


(2.9) 

( 2 . 10 ) 


Next, let iV > 2 be a fixed integer. We denote by 0 < j < TV, the zeros of the 
polynomial (1 — £ 2 )L' N (Q in increasing order. We recall (see [7, Chap. I]) that there exist 
positive weights pj, 0 < j < N, such that the following equality, called the Gauss-Lobatto 
formula, holds 


V$ e P 2 jv-i(A), f 1 *( 0 dC = 

' /_1 3=0 


( 2 . 11 ) 


Moreover, it follows from the identity (see [7, Chap. Ill]) 

£M0) 2 ft = (2 + ^)M, A , (2.12) 

j=0 


that the bilinear form: (u,v) —>• u (£j) v (£j)Pj is a scalar product on Pjv(A), since we 

have 

N 


Vv € Pjv(A), |M|o, a < £ v&fPj < 3||«||S, A . (2.13) 

l=o 


THE CONTINUOUS PROBLEM 

First variational formulation 

We define the subspace i7 1 (f};I v ) = {u € i7 1 (Q )\v = 0 on T'} of iL x (Q) and we 
introduce the space: 

M = H 1 (Q; T 2 ). (3.1) 

In the same way as in [1], we consider the following equivalent variational formulation 
of the problem (1.1)-(1.4): 

Find u in L 2 (Q) 2 and p in i7 1 (fl) such that p — ip belongs to M and that 

Vv e L 2 (fl) 2 , a(u,v) + b(v,p) = [ f(x).v(x)dx, (3.2) 

Jn 

\/q e Af, b{ u, q) =< U 0 , q > r i, (3.3) 

where the bilinear forms a and b are defined by 

V(v, w) <E (L 2 (ft) 2 ) 2 , a(v, w) — f v(x). w(x) dx, (3.4) 

Jn 

Vv e L 2 (Q) 2 , \/q G i7 1 (n), b(v,q)= [ v(x). Vg(x) dx. (3.5) 

Jn 


Theorem 3.1 Lei f be in L 2 (fil) 2 , Uq in H~ 1 / 2 (T 1 ) and (p in iL 1 / 2 (r 2 ), where H~ 1 / 2 ( r 1 ) 
and tfV 2 (r 2 ) are defined respectively in (2.3) and (2.2). Then problem (3.2), (3.3) has a 
unique solution satisfying 

|I u IIl 2 (o) 2 + Ibll/rqn) < C'dlfllz, 2 ^) + 11 £4) 11- 1 / 2 ,r 1 + IMIi/ 2 ,r 2 )- (3.6) 

Proof. First, let us define (p belonging to H l t 2 (F) such that <p ^2 = ip and ||(^||i/ 2 ,r < 
c|Mli/ 2 ,r 2 - To this end, we must extend ip to a function belonging to iL 1 // 2 (r). Let p be a 
function defined in [0,2] by 

p(t) m 1 — t, for 0 < t < 1 and p(t) = 0, for 1 < t < 2. 
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We define <p\r ZI by 

<f(ai + tm) = a 0 + (2 - t)ri) + fi{2 - t)ip(a n + (2 - t)T ln ) 

and (p\ Vlv by 

v( a o - tr IV ) = n(t)ip{ a 0 + trj) + - t)<p(a n + tr In ). 

Then we have ||^||i/ 2 ,r T C||V^lli/ 2 ,r 2 - Next, let <f> in H 1 (Q) such that <3?| r = (p. Finally, we 
obtain a function verifying 


$|r 2 — ip and ||^||ifi(o) < C1Mli/2,r 2 * (3.7) 

Second, let us extend Uo to a function belonging to H~ 1 / 2 ( r). We set C/oir 1 — Uo and 
Uo\r 2 = — \ < Uo, 1 >r x • Then, we have ||C/o|| —i/ 2 ,r < C'||{7o||-i/2,r 1 and < Uo, 1 >r= 0. 
Next, we define Neumann’s Problem: 


—= 0 

^ -u 

dn- U ° 


in Q 
on r 


and we set uq = Vt/>. Applying Proposition 1.2 of [11, page 14], we derive that ip belongs 
to H 1 ( f2), which implies that uq belongs to H( div; Q) with 


||uo||ff(div;n) — W h 1 ^) < C'||£?b||-i/ 2 f r> 


since div u G = Aip = 0. Finally, uq verifies 


u 0 .n| r i = Uq, divu 0 = 0 in Q and ||uo||jf(div ; n) < C||C7 0 ||_i/2,ri- (3.8) 

Now, let us split p as: p = + p with p in M and u as: u = Uq + u. Then, we can 

write the problem (3.2),(3.3) as 

Vv G L 2 (p) 2 , a(u, v) + b{y,p) = [ (f (x) — V$(x) — u 0 (x)). v(x) dx. (3.9) 

VqeM, b(u,q)= 0. (3.10) 

Since the right-hand side of (3.9) defines a continuous form on L 2 (Q ) 2 and since the prop¬ 
erties of continuity and ellipticity are obvious we have only to check the following inf-sup 
condition on the form b (see [11, pages 58,59]): 


inf sup ^ > (3 -<=>- \/q G M, 

9 eM vez, 2 (n) 2 |fy|U 2 (n) 2 


sup 

veL 2 (fi ) 2 


fr(v,g) 

MU 2 («) 2 


> P\\q\\HHn), 


(3.11) 


with a positive constant j3. This ” inf-sup condition” was introduced independently by 
Babuska [3] and Brezzi [9]. We can verify this condition by taking v = \7q . Indeed, we 
have 

b(y,q) ^ b(Vq,q) , , 

SU P FT- — ITvFT- = I'J|« 1 (S!) 

vei 2 (Si) 2 IMUw II vg||i2(n) 2 

and, since q^ = 0, using a generalization of Poincare inequality (see [11, Chap. I, page 
40]) yields ||(ji| L 2 (si) < V\q\w(n), which implies 


IklU 1 ^) - \](P) 2 + iJM-ffbii)- 

Thus, the ” inf-sup condition” is verified with the positive constant (3 = , — . Hence, 

yJW+ 1) 

applying Theorem 2.3 [7, pages 116,117], the theorem follows. <0 
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Second variational formulation 

We introduce the space: 

I = {vG Z/ 2 (Q) 2 ; div v G L 2 (Yl) and v . ri| r i 


0}. 


(3.12) 


In an analogous way as in [1], we consider the following equivalent variational formu¬ 
lation of the problem (1.1)-(1.4): 

Find u in X and p in L 2 (Q) such that u — Uq belongs to X, where Uo is the function 
previously constructed that verifies (3.8), and that 

Vv G X, a(u, v) + b*(v,p) = [ f(x). v(x) dx—< (/?, v. n > F 2 , (3.13) 

Jn 

\/q G L 2 (0), b*(u,q) = 0, (3-14) 

where the bilinear form b* is defined by 

Vv G H(div ; fi), \/q G L 2 (fl), b*(v,q) = — [ div v(x). q(x) dx. (3.15) 

Jn 

In the same way as previously, we split u as: u = Uq + u. Then, we can write the 
problem (3.13),(3.14) as 

Vv G X, a(u, v) + 6*(v,p) = [ (f(x) — u 0 (x)). v(x) dx— < p, v . n >p 2 . (3.16) 

Vq<EL 2 {Q), b*(u,q) = 0. (3.17) 

Since the right-hand side of (3.16) defines a continuous form on X and since the proper¬ 
ties of continuity and ellipticity are obvious, we have only to check the following inf-sup 
condition on the form b*: 

Vq € L 2 (Q), sup 6 (V ’ 9 ^ > /3*IMU 2 (si), (3-18) 

vGX 11 V 11 jtf(div ; f 2 ) 

\ r 

with a positive constant (3*. Let us note that qo = q — 7-77 / q(x) dx belongs to L%(Q) and, 

\u\ Jn 

owing to a classic result (see [11, Chap. I]), there exists vq in Hq(YI) 2 , such that 
div Vq = —qo and 

Then, we set 

v = v 0 + Vi with Vi(rr, y ) = 

We can verify that v belongs to X with 


l|vo||jyi(n ) 2 < c\\qo\\L 2 (n)- 
(-(T f n g(x) dx)x, 0), -1 < x, y < 1. 


div v = -q and ||v 1 1^(div ; o) < C||v|| H i(n)> < C"ll<?ll, 


since | q(x) dx| < ||?||.£»(n. Then, we have 


sup 

vex 


b*(v.q) 

MI//(div;S!) 


y(v, 9 ) 

v||fl-(div;!!) 


> 


^7 IMU a (n). 


Hence, we derive the inf-sup condition and we obtain the following result. 

Theorem 3.2 Let f be in L 2 (fi) 2 , Uq in /y _1 / 2 (|' 1 ) an d p in i/ 1 / 2 (r 2 ), where H 1//2 (r 1 ) 
and H l / 2 (Y 2 ) are defined respectively in (2.3) and (2.2). Then problem (3.13), (3.14) has 
a unique solution satisfying 


u||ff(div,n) + IHU 2 (n) < C(||f|| L2( a) + H^oll-i^r 1 + 1111/2,r 2 )- (3.19) 


22 





Bernard (2018) 


Regularity results 

When the data f is in H (div ; Q), taking the divergence of the first equation of the prob¬ 
lem (1.1)-(1.4) and owing to the other equations, we obtain a mixed problem of Dirichlet- 
Neumann for the Laplace operator: 


A p = div f 


P = 
dp 

dn 


V 

= f. n - U 0 


in O 
on r 2 

on r 1 . 


(3.20) 


We suppose that f is in if 1 )f!) 2 , tp is in ff 3 / 2 (r 2 ) and Uo is in H 1 / 2 ( r 1 ). In addition, we 
assume matching conditions at the vertices of r (see [7, Chap I]): 

[ - tTj) - (f .11 - u 0 )(a.j + f-r J+ i)| 2 ^ < +oo, J = LIII 

Jo dTj t 

[ \(f.n-U 0 )(a J -tTj)+ d ' iD (blj + tr j + i)| 2 ^ < + 00 , J = II,IV. (3.21) 
Jo dT J+ 1 t 


Theorem 3.3 For any data f in if 1 )!!) 2 , tp in ff 3 / 2 (r 2 ) and Uo in if 1 / 2 (F 1 ) , where 
ff 3 / 2 (r 2 ) and H 1 ' 2 ^ 1 ) are defined in (2.2) , verifying the matching conditions (3.21), 
the solution (u ,p) of the problem (l.l)-(l.f) belongs to iif 1 (Q) 2 x H 2 (Q). 

Proof. Owing to matching conditions (3.21), there exists po in H 2 (Q) such that p 0 |r 2 — y? 
and (^p^)ir 1 = f. n — Uq. Let us set p = p — po. The problem (3.20) is equivalent to the 
following problem: find p in H 1 ( fi; T 2 ) such that 


Vq G iL 1 (f2; T 2 ), a(Vp,'Vq) = [ (div f + Ap 0 )(xWx) dx. 

Jsi 

Since the boundary between T 1 and T 2 is the set of vertices of T, the regularity of the 
data implies that this homogeneous mixed problem of Dirichlet-Neumann for the Laplace 
operator has a solution p in H 2 (Cl) (see [12]). Hence, we derive the regularity of p and u. 


Remark 3.4 If (f. n— Uo) is Lipschitz-continuous onTj or belongs to iL 1 (Tj) 7 J = //, IV 
and if <p belongs to C 1,:L (r j) or to H 2 {Tj), J = I, III, the matching conditions (3.21) are 
equivalent to simpler conditions: 

= (f . n - U 0 )(a.j ), J — I, III and- — (a j) = (f. n - U 0 )(aj), J = 11, IV. 

dr j dr J+ , 

SPECTRAL DISCRETIZATION 

First spectral discretization 

We define the discrete scalar product by 

N N 

(u, v) N = “(&> QPiPj- (4-1) 

2=0 j= 0 

and we denote by In the Lagrange interpolation operator at the points (fi, £j), 0 < i, j < N 
in Pat(^). We set 


X w = P Ar (!!) 2 or Xn — (Pjv-i(A) ® Pat(A)) x (Pjv(A) <3 Pjv-i(A)). (4.2) 
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We assume that the data f belongs to C°(Q ) 2 and, for the sake of simplicity, that u and p 
satisfy homogeneous boundary conditions, that is to say, we set = 0, Uo = 0 in (1.1)-(1.4). 
From the variational formulation (3.2)-(3.3), we derive the following discrete problem: 
Find Ujv in Xn and p n in Pjv(^) D M, where M is defined by (3.1), such that 


Vvat £ Xn , (ujv, vn)n + &Ar(v.v,PAr) = (f, vn)n, 
VqN £ Pat(^) fl M, &jv(ujv, Qn) = 0, 


(4.3) 

(4.4) 


where the form 6 jv is defined by 

Vvat £ Pat(Q) 2 , VqN £ IPjv(^)i &jv(vat, gw) = (vjv, Xq ^ N - ( 4 . 5 ) 


We have a classical saddle point problem. We verify the inf-sup condition 

Vqn e Pjv(fi) n M, sup > 7 || 9jv || tfl(n)j (4.6) 

v N eX N ||^At||l 2 (fi ) 2 


where 7 is a positive constant independent from N, by taking vn = XqN- Hence, we derive 
the following theorem. 

Theorem 4.1 Let f be in C°(fl) 2 . Then problem (4-3), (4-4) h as a unique solution 
(u N ,p N ) satisfying 

ll u Jv|U 2 (S2) 2 + |lPJv||ifl(S!) < C ||Iwf||L 2 (SJ) 2 - (4.7) 

Next, we establish a theorem which implies the convergence of our discretization 
method. 


Theorem 4.2 Assume that the solution (u,p) of problem (4-3), (4-4) belongs to H S (Q) 2 x 
H S+1 (Q), s > 0, and the data f belongs to H a (Q) 2 , a > 1, where H s (Ll), for non-integral 
values of s, is defined in (2.1). Then, the following estimate holds 

||u—UAr||i 2( n )2 + ||p—pjvllff^n) < c (-V _s (||u||tf S (Q )2 + ||p||^+i(n)) + -^ _<T ||f||^(n) 2 ) • (4.8) 

Proof. From the abstract error estimate for the approximation of saddle-point problems 
(see [7, Chap. IV]), we derive the following estimate: 


+ inf 

vweXjv 


u — ||z, 2 (o) 2 + lb — Pivlljprn) < c ( inf ll u — w ivlU 2 (Si) 2 

\wn€Vn 

In vjv (x) . zjy (x) dx - ( v N , z N ) N 
ll z Jv|| L 2 (S!) 2 

. ( /I, I, . b(z N ,q N ) - b N (z N ,q N ). 

iSLnJWv - + _ SU P - 


- Viv||L 2 ([i ) 2 + SUp 
zn£Xn 


9NGlPiv(fi)nM 


zjveAw ll z ivlU 2 (S2) 

fsi f(x). Zjv(x) dx - (f, zjv)jv 
+ SU P -n—n- 

W <EX n ll z JV||i 2 (S2) 2 


where Vjv is defined by 


(4.9) 


Vjv = {w w € Xn', VqN 6 Pjv(fl) n M, b N (-w N ,q N ) = 0}. 

Moreover, we recall (see [11, Chap. II, (1.16)]) that 

Q 

inf ||u - wjv||l 2 (! 2 ) 2 < - inf (||u - vjv|U 2 (Si) 2 ). 

w w evy 7 v N ex N 

Hence, we derive that, for all vjv_i and fjv-i in Pjv-i(H ) 2 and all q\- i € II Vy i(Q) n M, 
ll u — u v|U 2 ([i) 2 4- ||p — pjv||rfi(n) 

< c(||u — vjv_i||L 2 (n ) 2 + ||p — gjv-illffUn) + ||f — fjv-i|U 2 (si ) 2 + ||f — 2ijvf||L 2 (!j) 2 )- (4.10) 
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Then, we choose vjv-i ~ Iljv-iu (resp. fjv-i = n^-if), that is to say the orthogonal 
projection of u (resp. f) on Pjv_i(Q) 2 in L 2 (Q) 2 and qN-i — n]v- 2 iP? where is the 

orthogonal projection of p on Pn_i(£ 2) fl M in H 1 (H). It remains to prove the estimate, 
for any m > 1, 

Vp e |b-n^ lP ||^ (!3) < C Ar 1_m |b|| ff m ( n). (4.11) 

On the one hand, this result is obvious for m = 1. On the other hand, for ra > 2, we have 
(see [7, Chap. Ill]), 

lb - njv^iPlb'w - ID inf Wv ~ 

7’jv_ielFV_i(s2)nM 

< lb _2 JV-lP|k 1 (!)) < C A rl_m |blU m (!!)- 

Then, an interpolation argument (see [11, TH 1.4, page 6]) gives (4.11). Finally, the re¬ 
sult follows from (4.10), (4.11) and the classic estimate for the orthogonal projection on 
Pjv_i(fi) in L 2 (P). 0 

Remark 4.3 With the choice Xn = Pn{Q) 2 , problem (4-3, (4-4) can be interpreted as a 
collocation scheme. Indeed, by integrating by parts in the discrete bilinear form bjsr with 
respect to one of the two variables for each of the two terms of b^ (this process being 
allowed by the precision of the quadrature rule), and choosing as test functions the Lagrange 
polynomials associated with the grid points of Sjv, it is easily seen that (4-3), (4-4) 
equivalent to the set of equations for u n in Pn(Q) 2 and pn in Pn(LI) 0 M: 

Ujv(x) + Vpat(x) = f(x), Vx G Sjv, 

divujv(x) = 0, Vx G S/v H 12, 

2 

divu w (x) = (ujy . n)(x), Vx€=jvnr. 

Second spectral discretization 

In order to improve the approximation of the condition divu = 0, we can try to 
decrease the dimension of the space X^. So, we choose 

X N = PWfi) 2 . (4.12) 

We note that, in this case the forms b (.,.) and 6 jv(-> •) are equal on Xn x Pjv(^)- It does 
not appear spurious modes for the pressure, as we can see in the following lemma. 

Lemma 4.4 Let Zn be the space 

Zn = {qN € Pn(£1) L\ M; Vv^r G Pjv-i(^) 2 , 6 (vjv, #at) = 0}. 

Then Zn = {0}. 

Proof. Let qjy be in Z^. Since j s a polynomial of degree < N — 1 with respect to x, 

ox 

which is orthogonal to Pjv-i(^), we can write: 

q N (x,y) = a N (y) + p N (x)L N {y), 

with ckjv and (3n in Pjv(£ 2). In the same way with y in place of x, we have: 

q N {x , y) = 7 n(x) + S N (y)L N (x), 
with and Sn in Pn(LI). Hence, we derive 

q N (x, y) = A + /iL N (x)L N (y ), 
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where A and p are real numbers. But, since L^( 1) = 0, the condition: 

Vy € [- 1 , 1 ], gjv(l, 2 /) = 0 

encienciaje 

ingenierias 

implies A = p = 0. 

We have to study the following discrete problem: 

Find Ujv in Pat-i ^) 2 and pn in Pat(^) D M such that 

Vv^ € Ptv-i(^) 2 , (u^,vat)at + ^(vat,Pat) = (f, vat)jv? (4-13) 

€ Pat(^) H M, b(uN,qN) = 0. (4-14) 


aci 


For the inf-sup condition, the choice v/v = 'Vqn is no longer available, since v^r must 
be in Pjv-i(n) 2 . In fact, in the next lemma, we find a inf-sup constant depending on N. 

Lemma 4.5 There exists a constant c > 0 independent on N such that 

Vq N e ZPjv(O) n H 1 ^; r 2 ), sup .f (v ^ ,gJr) > cATliav-llHHn)- (4.15) 

v N eiPN-i(fi ) 2 II v at||z, 2 ( 0) 2 

On the one hand, we note that 

N 

<&(. x >v) = a n (x)L n (y), 

n=0 

where a^ is a polynomial of degree < N — 1. Then, we have 

^r( x < y ) = n AT-i(^)(a:, 2/) + a' N {x)L N {y), 
which implies, using (2.7) and the inverse inequality (2.10), 

||-J^||L 2 (fi) < ||n iV -i(-^)||L 2 (n) + ciV2||a J v||o,A- (4-18) 

On the other hand, in view of (2.8), we can write 
da* N 

n N -i(^)((a:, 2/) = ^ nAr_i(a n L^)(x, y) = (27V - l)ajv(a;)Lj\r_i(y) + rjv(a;, y), 

Oy n=0 

where rx(x,y) is a polynomial of Pjy(fi) of degree < N — 1 with respect to y. The 
orthogonality properties imply 


||njv-i(-^)||L 2 (fi) > 2y TV - -||ajv||o,A- 
By combining this inequality with (4.18), we obtain 

II^IU 2 (!i) ^ ll n JV-l(-|^)llL2 (!i) +cA r ||n J v_l(-^)|| L 2 (n) . 


(4.19) 


9Qn 


In the same way, we have the analogous inequality for ——. Thus, we obtain 

ay 

II ^QnWl 2 ^) 2 < cJV||n w _i(V3^)||i?(n)». 

Next, the equality (4.16) yields 


(4.20) 


IIVqrjv[|i=(n) < 1^17^11x2(5,, + |a J v|||V(LAr(^)LAr(j/))||L 2 (!i): 
which implies, owing to (2.7) and (2.9), 


26 
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I|V?rv|U 2 (!i) < l|Vg^-||L 2 (ti) + 2y 9j\r + 1 ^ l ajv l' 

It remains to estimate |ajv|. First, we have, owing to (4.16), 

N 

q N (x, y) = Y1 a n {x)L n (y) + a N L N (x)L N (y). 

71=0 

But, V?/ € [—1,1], 2 /) — 0, which implies 


a n (T) = 0, n = 1,..., N — 1 and ajv(l) + ojv = 0. 


(4.21) 


N—l 

If we set a^(x) = ^ a k L k (x), we derive 
k =1 


N-l 

a N = -^2 a k , 
k =o 

and, therefore, thanks to a discrete Cauchy-Schwarz inequality 

M < (E ^^E^ + h )) 1 <^n\M\ 0A . 

k=0 K 2 fc =0 Z Z 

Then, in view of (4.19), we obtain 

\&n\ < 'C l|IIjv_i(Vg^)|| L 2(^). 

Finally, (4.20), (4.21) and (4.22) yield 


(4.22) 


||V<7A/-||z,2(n) 2 < clV||n 7 v_i(V^)|U 2 ( fi ) 

which, in view of (4.17), ends the proof. <0 

The bilinear form .)n and b(.,.) satisfy Brezzi’s conditions with respect to Pjv-i(f ^) 2 
and Pjv(Q) n M (the bilinear form (.,.)n is continuous on Pjv-i(f ^) 2 and elliptic on 
Pn-i(^), the bilinear form &(.,.) is continuous on Pjv-i(f^) x (Pjv(Q) fl M) and verifies 
the “ inf-sup condition”), see [7, Theorem 2.3, pages 116,117], whence the theorem. 

Theorem 4.6 Let f be in (7 0 (n) 2 . Then problem (4-13), (4-14) has a unique solution 
(u N ,p N ) satisfying 

|| u iv|U 2 (ft ) 2 + N -1 ||pAr||/fi(n) < C \\Tn^\\l 2 {Q) 2 - (4.23) 

We establish the convergence of this second discretization in the following theorem. 

Theorem 4.7 Assume that the solution (u.p) of problem (4-13), (4-14) belongs to H S (Q) 2 
xH s+1 (Q), s > 0, and the data f belongs to H a (Q) 2 , a > 1. Then, the following estimate 
holds 

Il u — u Ar||z, 2 (ft ) 2 + N -1 ||p — < c (Ar -S (||u||tfs ( fi )2 + ||p||ff«+i(n)) + N -<T ||f||^(n) 2 ) • 

(4.24) 

Proof. The abstract error estimate, analogous to (4.39) but with much simplification 
because of the exactness of the quadrature formulas, yields 

II u - ujv||l2(!j)2 + lV _1 |b - < C ( inf ||u - WjvIIms !) 2 

\w n€Vn 
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+ , D inf , n ., ll u ~ v iv||i, 2 (Si) 2 + inf 1/f l|p-9Jv||ffi(£!) + ||f-^Jvf|U 2 (n) 2 , 

vjvGlPjv-i(r2) 2 g;\rGlPjv(fi)nM J 

encfencfale where Vn is now the space 
ingenierfas 

Vn = {w N € Fat-i(^) 2 ; 6 Pjv(^) D M, 6(wat, <?at) = 0}. 

It remains to estimate the term inf ||u — Wjv||z, 2 m) 2 . Since u is such that divu = 0 and 

wjveVjv 

u. nipi = 0, there exist a unique ip in H 1 (H) (see [11, Chap. I and 4]) such that 
u = curl^> and ip = 0 on T 1 . 

Moreover, if u belongs to H S (Q) 2 , we have ||'0||/a rs + 1 (n) < c ll u lli? s (n) 2 - Let us define the 
operator 7f^ (see [7, Chap. II]) on if 1 (A) by 

V<p 6 H\ A), (*»(C) = (4V)(C) + ¥>(-1)1^ + ¥>(1)1^, (4.25) 

where the function <p stands for 

m=v(o - v>(-i)b=^ - aV-V- 

Note that the definition of is available, because (p belongs to i7g(A), and that 7rJ) <p and 
tp coincide in —1 and 1. In [8, Section 7], the following estimate is proven, for all r > 1 
and all function <p in H r ( A): 

\<P - 7Tv^|l,A + N\\<P - 7f^||o,A < C A^ 1_7, ||^|| r)A . (4.26) 

Assuming s > 1, we set 

iijv-i(u) = curl o 7f^ 1 V')- 

Since ( 7 fj,*, o i r i = j/'ir 1 — 0, we can verify that ll\- ](u) belongs to Vn and, in 

view of (4.26), that 

W-t^n-x 0 7i'v s) iV’|if 1 (si) ^ ciV~ s ||V’ll/P+i(Si)- 
Finally, we derive, for s > 1, 

w ^nf^ ||u — w„|| iW < ||u — i?Ar-i(u)||L 2 (ft) 2 < cJV-Hull^. 

Since we have inf ||u — Wjvll^afma < c||u||jr, 2 / n ), an interpolation argument gives the 

w n€V n 

result of approximation in Vn for any s > 0 and the estimate of the theorem follows. <() 



Third spectral discretization 

The third discretization comes from the variational formulation (3.13), (3.14). We 
define the space Xn by 

Xn = Pj\r (^) 2 n X = {vjv £ P;v(^) 2 ; vn ■ n\ri = 0 }. 

Let Mn be a subspace of Ptv(^) that we shall set later. Then, we consider the following 
discrete problem: 

Find ujv in Xn and pn in Mn such that 

Vvat € Xn, (ujv, vn)n + b* N (vN,PN) = (f, vn)n, ( 4 . 27 ) 

WqN £ Mn, b* N (\iN,qN) = 0 , ( 4 . 28 ) 


28 
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where the form b* N is define by 

Vvjv € Pjv(n) 2 , Vgjv € Pjv(O), b* N (v N ,q N ) = -(divv N ,q N ) N . (4.29) 

In order to choose Mn, we begin to identify the spurious modes for the pressure. These 
spurious modes for the pressure are derived by elimination from those of classic Stokes 
problem (see [7, Chap. IV]). In particular, we can verify: Vv/v G b* N {vpj, q^) = 0, for 
qN(x : y) = Ln(x) or Ljv(a;)Ljv(?/). We obtain the following lemma. 

Lemma 4.8 Let Zn be the space 


Z* N = {qN G Pn(Q)', Vvjv G Vat, b* N (vN,qN) = 0}. 

Then is spanned by (Ln{x), Ln{x)Ln (y))- 

Finally, let Mn stand for the orthogonal complement of Zn for the scalar product in 
L 2 (Q) or for the scalar product (., .)n, owing to (2.11). The inf-sup condition is given in 
the next lemma. 

Lemma 4.9 There exists a constant c > 0 independent from N such that 

\fqN G M/v, sup n( n,Qn) > c \\q N \\ L 2 (n). (4.30) 

vatGXjv ||Wjv||.ff(div;n) 

Proof. Any function qN in Mn has the expansion 
N-1N-1 

qNfay) = X X qm,nL m {x)L n (y) 

m=0n=0 

N- 1 N-l 

+ q™,NLm(x){L N (y) - L N _ 2 (y)) + Y1 <lN,n{ L N(x) ~ ijV-2 (x) )L„(y). 

771=0 77=1 


With the convention L_i = 0, we choose wn = (wn: zn) with 


2m + 1 


, , L m+1 (x)-L m _ x(x) 

w N (x,y) = - Xs N 9™."- , , - L, l {y) 

771=0 77=0 

L m +i(x) — L m _i(x) , . . / \\ 

— X 9 m,JV-, 1- (LN{y) - L N -2(y)) 

777=0 


2m + 1 


and 


7 x ^ ^ ,7 .L n+1 (y) - L^y) 

Zn{X, y) / y / y qm,nLm\X) 

777=0 77=777+1 

AT—1 

X QN,n{LN(x) — L N _ 2 (x))~ 


2n + 1 
L n+ i(y) - L n _i(y) 


2n + 1 

Then, in view of (2.8), we have 

divwjv = —qN and w n G Xn 

since w n ■ nipi = ^jvir 1 — 0. As in [6, Chap. IV], we prove 

/N-1N-1 i 


IMIi* ( n) >c EE« 


'"(m+|)(n + |) 

i N-l i N-l 

+ ^,AT^T + X VN,r. 

iV ' o m—n iib-r o „_i 


n + 


t) 


(4.31) 


(4.32) 


(4.33) 


(4.34) 
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Hence, we derive, in the same way as in [8, Section 24] 

ll-^IU^S!) + ll-T^IU^n) < c||gjv|U=(n)- (4.35) 

Next, setting w^(x : y) = wn(x, y ) — Qo,oLi(x) — qo,NLi(x)(LN(y) — -^- 2 ( 2 /)), we note that 
tujr(±l, y) = 0 for — 1 < y < 1. Then, the Poincare-Friedrichs inequality, applied with 
respect to x or y, yields 


II“'5 V || L 2(S2) < c||-^||l2(sj) 
Hence, owing to (4.35) and the estimate 


and 11zw11z, 2 ( fi ) < c||-^|| L 2 (ft). 


\ 


<&,0 


Qo,n 

n T 


T - c I|9n||l2(!2), 


(4.36) 


which is derived from (4.34), we obtain 

||Ww||i 2 (n) + ||zjv||.W(n) < c ||<?jv||i 2 (! 2 ). 

Finally (4.35) and (4.36) imply 

||wjv||iT(div;!i) < c||gjv||L 2 (n) 

and, in view of (4.33), the choice vjv = wjv in ^ v (vv, qx) is available and gives the inf-sup 
condition (4.30). <0 


From the previous lemma, we derive the following theorem. 

Theorem 4.10 Let f be in C°(fl) 2 . Then problem (4-27), (4-28) has a unique solution 
(ujv, Pn) satisfying 

ll u iv||ff(div ; s!) + ||pw|U 2 (o) < OHXjvflli.^.. (4.37) 

Sketch of the proof. From(4.27), we derive (ujv,ujv)jv = (T.vf. u.vj.v- Then, Owing to 
divujv = 0 and (2.13), we derive |uv||^(di V ; 0 ) < 3||X\-f | f, 2 (n)■ Next, the inf-sup condition 
(4.30) and (4.27) imply 


lbjv||r,2(!2) < 


C vjvCXjv |Mlff(div;S2) 


1 ( 2 )vrf, vjv)at — (ujv, v N )jy 

- sup -jr n - 

^ vj yEXpi ||^Ar||i/(div;f2) 


Hence, in view of (2.13), we obtain 

18 

lbw|U 2 (Si) < — ||2)vf||L 3 (n), 


which ends the proof. <0 

Next, in the same way as in Theorem 4.2, we prove an optimal error estimate. 


Theorem 4.11 Assume that the solution (u,p) of problem (4.27), (4-28) belongs to H S (Q) 2 
xH s (Q), s > 0, and the data f belongs to H a (Q) 2 , a > 1. Then, the following estimate 
holds 

Il u_u Ar|U 2 (ft) 2 + \\p~P n\\l 2 (D) < c (-^ -S (ll u ll/r s (n) 2 + ||p||tf s (n)) + -N _CT ||f||/z ff (ft) 2 ) • (4.38) 

Sketch of the proof. Again, from the abstract error estimate for the approximation of 
saddle-point problems (see [7, Chap. IV]), we derive the following estimate: 
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||u — Ujv|| L 2 (si )2 + ||p-piv|U 2 S!) < c ( inf ||u - WjvHz^S!) 2 

\w n€Vn 

J n vjy(x) .Zjy(x)rfx- (vjV,Z]y)jV 

ll z JvlU 2 (£l) 2 

b*(z N ,q N ) - by(zj y,q N ) 

II z ;v||l 2 (!2) 2 


+ i°t (||u - vjv|U=(m 2 + sup 
v N £X N zn£X n 

+ inf (||p-9jvIU 2 (sj)+ sup 


+ sup 

z n€Xn 


In f (x). zjy(x) clx - (f, zjv)jy 

ll z w|U 2 (!2) 2 


(4.39) 


where Vjv is defined by 


Vn — {wjv € \/qx € Mjv, b* N (\VN,qN) — 0}. 

Moreover, we still have (see [11, CH. II, (1.16)]) 

Q 

inf Hu-WjvIU 2 ™ 2 < - inf (||u - vjvll^rm 2 )- 
w N ev N 7 v N ex N 

We end the proof in the same way as in Theorem 4.2. 4} 

Remark 4.12 As for the first discretization, problem (4-27, (4-28) can be interpreted as 
a collocation scheme. In the same way, we prove that (4-27), (4-28) is equivalent to the 
set of equations for in Xn and p^ in 


Ujv(x) + Vpiv(x) = f(x), VxeSjvflQ, 

N(N + 1) ( (UJY ' + % f = (f ' n)(x) ’ Vx € En n r2 ’ 

divujv(x) = 0, Vx G Sjv- 


Therefore, the discrete solution u^r is exactly divergence-free, which is important for 
some applications. 


NUMERICAL RESULTS 


The convergence of the methods corresponding to the first and third discretizations 
were tested in a problem of the type (1.1)-(1.4), with homogeneous boundary conditions. 
Precisely, we tested the convergence of these methods to the exact solution u(x, y ) = 
(ttx 2 cos(7r?/), — 2xsm(ny) and p(x,y) = t/sin(7nr), which means that we stu¬ 
died the convergence of these methods for Problem (1.1)-(1.4), with 

f(x, y) = (ttx 2 cos(7n/) Try cos(ttx) ,—2x sin(iry) + sin(7rx)) 

and homogeneous boundary conditions. In addition, we tested the convergence to 0 of the 
divergence for both methods. 

We shall use the Lagrange polynomials. We denote l r the Lagrange polynomial asso¬ 
ciated to the Gauss-Lobatto point 0 < r < N, the expression of which is 

n (*-&) 

lr(x) = ^-• (5.1) 

n (fr-fj) 

J=0, jj^r 

The derivative l' r verifies the following equalities 
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_n (fm-W 1 

Vr = 0,..., JV, Vm = 0,...,JV, r^m, £(£„,) = —|r"-(^—^r), (5.2) 

n (Cr-C#) u 4r 

j=0,j^ T 


Vr = 0 ,...,N,U(r)= £ (5.3) 

j=0,j^r 


Uzawa's algorithm 

Problem (4.3), (4.4), Problem (4.13), (4.14) and Problem (4.27), (4.28) are equivalent 
to a linear system of the type : 

f MU + DP = MF, 

[d t u = o ■ (54) 

The unkowns are the vectors U and P which represent respectively the velocity and the 
pressure values on a given grid points. The data / is representated by the vector F on 
the same grid points. The diagonal matrix M is the weight matrix, while the matrix D is 
associated to the form b n for the first and second spectral discretizations and to the form 
for the third spectral discretization and D T is the transposed matrix of D. 

Uzawa’s algorithm consists in rewriting the first equation of system (5.4) as: 
U = F — M~ X DP and substituting in the second equation. We obtain a new equation for 
the pressure P : 

(D t M~ 1 D)P = D t F. (5.5) 

Next, we solve this symetric system either directly if the matrix D T M~ 1 D is invertible or 
by diagonalizing the matrix D T M~ 1 D if not, because the spurious modes correspond to 
eigenvalues of the matrix equal to zero. Next, we compute the velocity via the formula 

U = F - M~ l DP 

and its divergence by multiplying U on the left by the matrix D T . Finally, we test the 
convergence to 0 of its divergence. 


Implementation of the rst discretization 

We take (lj(x)lk(y),0), (0 ,lj(x)lk(y)), 0 < i,j < TV as a basis of Xn and 
l r (x)l s (y)i 1 < r < N — l,0<s<TVasa basis of Pn{&) FI M. The unknowns are 
the velocity un = {un^Un) and the pressure p n. For j = 1,2, for 0 < r, s < N, we 
denote = u J N (£ r ,£ s ), F^X = fj(£ r ,£ s ), where f = (/ 1} / 2 ) represents the data, and, for 
1 < r < N — 1, for 0 < s < TV, we denote P^ s = pjv(?r, Cs)- So, we have 

u J N (x, y) = U>fl r (x)l s (y), j = 1,2 and p N (x,y) = Y1 Tj P r,s l r(x)h(y)- (5.6) 


Let us define the (N + l) 2 x ( N + l) 2 diagonal matrix M = (^(j,fc),( 7 -,s))o<j,A:,r,s<iv with 




0 if (r, s) (J, k) 
PjPk if (r,s) = ( j,k ), 


(5.7) 


and the 2(JV + l) 2 X 2(N + l) 2 diagonal matrix 


M 


M 0 \ 

0 M )' 


32 


By setting \jsi(x,y) = (lj(x)lk(y), 0), for 0 < j, k < N and, next, vw(x,y) = (0 ,lj(x)lk(y) 
in (4.3), (4.4), we obtain the matrix system (5.4) where the column matrices U, P and F 
are defined by 
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U= ( U j£, ) , P= ( P r N s ) and F= ( jjjj j , 0 <j,k,s < N, 1 <r < N- 1, 
and, with (., ,)jv defined by (4.1), the 2 (N + l) 2 X ( N 2 — 1) matrix D by 

D= (y d 2 ’ n ) with DlN = ° 2 ' N - {(h( x %{y)M x )i's{y)) N )i 

where (j, k ) represents the row index, (r, s) the column index, 0 < z, j, s < TV, 
1 < r < TV — 1, for the matrices .D 1,Ar and Z) 2,Ar . Note that U and F are 2 (TV + l) 2 x 1 
matrices and P is a (TV 2 — 1) x 1 matrix. 


Proposition 5.1 The (TV 2 — 1) x (TV 2 — 1) square matrix D T M 1 D is invertible. 


Proof. We just have to prove that the rank of the matrix D is TV 2 — 1. Let us assume 
that the rank of the matrix D is strictly smaller than N 2 — 1. Then, there exist a sequence 
of real number (q r , s ), 1 < r < TV — 1, 0 < s < TV, where all the real numbers q r>s are not 
equal to zero, such that 

N-l N 


E E^Ay, = o, 


N-l N 

where D rs are the column vectors of the matrix D. Setting q = Qr,slr{^)h{y) , the 

r=l s=0 

previous equality is equivalent to 


Vvat G Xn, b(vN,q) = 0, 

which is in contradiction with the property that there is no spurious mode. <0* 

We have to compute the matrix D T M~ 1 D = (b(t, u ),{ r ,s)) ■> 1 < r, t < N— 1,0 < s,u < N. 
Owing to the previous expression of the matrices D and M, we obtain 

N 1 

b(t,u),(r,s) ^ y (/j7i(x)/fc(?/),/j(x)/ w (2/))jv(^m(‘^')^fc(2/)?^ r (^')^s(2/))jV 

m,k=0 PmPk 

N X 

+ E -(■ lm( x )lk(y),k( x )l' u (y))N(lm( x )lk(y),lr(x)l' a (y))N■ 

m,k= 0 

Next, we change the numbering for the matrix D T M~ X D. Let us define the mapping 
T> by 

V(r, s), 1 < r < N — 1, 0 < s < TV, cp(r , s) = 1 + (r — 1)(TV + 1) + s. (5.8) 

Note that is a one to one mapping from {1,..., TV — 1} x {0,..., N} to {1,..., TV 2 — 1} 
and we note 

VI < i < N 2 - 1, <^ _1 (z) = (-01 (i),V>2(i))- (5.9) 

Note that ^(z) — 1 and ^O) are respectively the quotient and the remainder of the 
euclidian division of z — 1 by TV + 1. Then, we can denote 

D M. D {'Q'i,j')i<i,j<.N 2 —i with cijj b(t,u),( r ,s )> (5.10) 

where (i,w) = (^i(z), ^(i)) and (r,s) = (^iO), ^O))- 

Computing the elements of the matrix D T M~ 1 D yields 

1) if ^i(z) ± M) and ± foti) 

,j 0 

2) if Vh(0 ^ VhO) and ^ 2 (i) = ^ 2 O') 

a i,j = Pip2(i) X) Pm0 1 (i)(Cm)0 1 (,)(Cm) 

m=0 

3) if Vi(*) = V'i(i) and tp 2 (i) ± M) 

N 

a i,j = Pi’l (0 S^/ ? m^ 2 (i)(^m)^ 2 (j)(Cm,) 
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4) if ipi(i) = (j ) and ^> 2 (0 — ^2 (j ) , that is to say i=j 

AT AT 

a i,i — 2 (i) Pm{J j xh i (j)) T Pip^i) X^ Pmil^U)) • 

771=0 771=0 

Thus, most of elements of the matrix D T M~ 1 D are equal to zero. 
Next, we determine the column matrix D T F = (c* } i) by 


t — N 1 , C i( i Pip 2 (i) -^771,1/12(7) PV’lf*) Pm^ 2 (i)^ m ')-^ , ipi(i), , rn- 

771=0 771=0 

Hence, owing to (5.2) and (5.3), we compute the list [//;(£ m ), 0 < k,m < JV], which 
allows us to determine the elements of the matrices D T M~ X D and D T F. Since the matrix 
D t M~ 1 D is invertible, the equation ( D T M~ 1 D)X = D T F has a unique solution X. Then 
we derive easily the column matrix P = (P^ 8 ) such that P^ s = X^( r>s ), 1 < r < N — 1, 
0 < s < N and, next, the column matrix U, thanks to the relation U = F — M~ l DP. 
Setting Pq ± = P$ t = 0 for 0 < t < N in accordance with the boundary conditions, we 
obtain 

- E u 2 m N k = - E m k )p^ t . 

11=1 t =0 


r / 1 - N . 

u m,k ■ 


: F 1,N - 

m.k 


Finally, we derive 


(div U^, div U N ) N = pipj{J2 ( U mj C( 6 ) + Cfe ))) 2 

*J=0 771=0 


(u - u N , u - u N ) N = y 5,) - I#*) 1 + (uate, 4) - Ulff), 

i,j =0 

(p-pn,p-pn)n = E E w(p(&.&) - ^j) 2 - 

1=1 j =0 


We give the values of (div u^v, div Ujv)jv, (p — Pat, p — Pn)n and (u — U/v, u — ujy)^- for 
N between 4 and 21. 


N 

4 

5 

6 

7 

8 

9 

T 

(div Ujv, div Ujv)jv 

6,85 

0,27 

1,12 

0,0155 

0,080 

5, 16 .IO - 4 

1 

*<3 

7 

T3 

0,06 

0,02 

2, 7.10 -3 

7, 75.10- 4 

7,59.10- 5 

1,74.10-® 

(u— ujv,u— ujv)|r 

1,07 

0,045 

0,010 

1,56.lO -3 

4,65.10- 3 

3,52.10-® 

N 

10 

11 

12 

13 

14 

15 

1 

(div ujv, div Ujv) jy 

3,14.10 3 

1,11.10 5 

7, 86 .IO - 5 

1, 70.10- 7 

1,36.10 6 

1.92.10- 9 

{jP-Pn,P-Pn)s 

1,44.10~ 6 

2, 78.10~ 7 

1,96.10 -8 

3,27.10- 9 

2,03.10- 10 

2,97.10- 11 

1 

(u-u N ,u-u w )^ 

1,33.10 -4 

5, 6 .IO - 7 

2,56.10-® 

6,62.10 9 

3,54.1Q- 8 

6,03.10-“ 

N 

16 

17 

18 

19 

20 

21 

(div u N , div ujv) h 

1, 73.lO -8 

1,68.10 11 

1,69.10- 10 

1,02.10- 12 

3, 79.10~ 12 

1,91.10-“ 

(j>-Pn,P~Pn)n 

1,64.lO -12 

2,15.10~ 13 

1,06.10- 14 

5,50.10~ 15 

1.09.10- 14 

1.06.10- 14 

1 

(u-u w ,u-ujv)a, 

3, 70.10~ 10 

4,37.10 -13 

3,03.IQ -12 

3,03.IQ -14 

9,03.10- 14 

2,69. IQ- 13 


Table 4.1 


We shall comment these results on comparing them with these ones of the third dis¬ 
cretization. 
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Implementation of the first discretization 


We only sketch the method because the first and third discretizations are rather similar, 
but we shall point out the differences between both discretizations. 

First, the space Xn is not the same, because boundary conditions are taken in account. 
Thus, we take (lj(x)lk(y), 0), (0, l s (x)l r (y)), 0 < i,j,s < N, 1 < r < N — lasa basis of 
Xn- 

Second, we shall compute the pressure p^ as the orthogonal projection on Mn of an 
element of P^(^) and we take l r (x)l s (y), 0 < r, s < N as a basis of For 0 < r, s < 

N, we denote = Ux(£ r ,£ s ), = p.v(£r,fs) and for 0 < r < N, 1 < s < N — 1, we 

denote = u%(^ r ,^ s ). Let us define the matrix Mi = M, where M is given by (5.7), 
the (N 2 — 1 ) x (N 2 — 1 ) matrix M 2 = { m (j,k),(r,s)) with 


^(j,fc),( r,s) 


0 if (r, s) ^ (j, k) 
PjPk if (r,s) = (j,k), 


0 < j,r < N, 1 < k, s < N — 1 


and the 2N(N + 1) X 2N(N + 1) diagonal matrix 


M, = 


Mr 0 \ 

0 M 2 )■ 


In the same way as the first discretization, we derive the following matrix system 


M.U + D,P = M t F, 

D^U = 0 


(5.11) 


where the column matrices U, P and F are defined by 

-U ,N 


U - 


U]f 

ut N 


F, 


t,u 


, P = (P^) and F = ( j , 0 < j,k,r,s } t < N, 1 < u < N — 1 


and the 2N(N + 1) X (IV + l ) 2 matrix D„ by 


D, = 


n with D\ ,n = (-(l'Ax)l k {y),l T (x)l s (y)) N ), D 2 t ' N = {-(l t (x)l' u (y) : l T (x)l s (y)) N ) 


We have to compute the matrix 0 < r, s,t,u < N . In the same 

way as the first discretization, we obtain 

N 1 
m,k =0 PmPk 

N N-l i 
m=0 k=l PmPk 

Next, as in the first discretization, we change the numbering for the matrix Dj M^D*. 
Let us define the mapping ip* by 


V(r, s), 0 < r, s < N, <£>*(r, s) = r(N + 1) + s + 1. (5-12) 

Note that ip* is a one to one mapping from {0,..., N} 2 to {1,..., (N + l) 2 } and we note 

VI < i < N 2 - 1, <p~\i) = W(i)i 2 *W)- ( 5 - 13 ) 

Note that (i) and are respectively the quotient and the remainder of the euclidian 
division of i — 1 by AT + 1. Then, we can denote 

= Kj)i<i,j<(iv+i ) 2 with aX = b* MM , (5.14) 
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where (t,u) = W’J(i), V’sKO) and (r,s) = (ipl(j), ipZU))- 

Computing the elements a\. } of the matrix I)JM~ l D t yields 

1) if ipl(i) =£ and V>J(i) f Va( 3 ) 

<,i = 0 

2) if ± and y>|(*) = 

n 1 ( 

a i,j = P$l<f)Ptl>t<j)P<l>Z(i) S 

?71=0 

3 ) if y$(i) = V’iO') and V’K*) ^ V’Ki) 

771=1 

4) if y£(i) = V’lO") and V’2*(i) = V’l(i). that is to say i=j 

AT AT—1 

a i,i = (Pi/>l(i)) PipZ(i)) (*)) (^ 2 (*))) ' 

zero. 

Next, we determine the column matrix DjF = (c* x ) by 


® — (-^ + 1) ) c i,l — P'01 (i)/ ? V , 2(*) ( L&iW)^rn',^(i) ~b ^fc(^2(*))^i(*)i«i)‘ 

771=0 771=1 

Now, we deal with the main difference between both discretizations. In the third 

discretization, the matrix Dj M~ l D* is not invertible and the computation of the pressure 

AT 

is more complicated. First, let = ^2 Qt,uh( x )lu(y) be a spurious mode, that is to say 

t,u=0 

an element of the space Zjy defined in Lemma 4.8. In the same way as in the proof of 
Proposition 5.1, considering that the matrices DjM~ l D^ and D * have the same rank, we 
obtain the following equivalences 

q N eZ* N ^=> D*Q = 0 (DfM-'DJQ = 0, 

where Q is the column vector (qt, u )o<t,u<N- We can consider the matrix D^M~ X D^ as the 
matrix of a linear mapping / from the vector space Pjv(^) into itself equipped with the 
basis B = (h( x )lj(y))o<i,j<N • Therefore, Z* N is the eigenspace associated to the eigenvalue 
equal to 0 of the linear mapping / or, equivalently, of its matrix DjM^D*. Note that 
this basis is orthonormal for the scalar product (•, -)n- We can diagonalize the positive 
symmetric matrix DjM^D* and, thus, there exist a diagonal matrix A and an orthogonal 
matrix R such that D* Mp 1 D* = RAR -1 and such that the diagonal elements of A, that is 
to say (A i)i )i< i <(jv + i )2 are in increasing order. Note that the matrix A is the matrix of / in 
a basis B' = (hi)i<i<(N+i) 2 °f Pat(^) 5 which is orthonormal for the scalar product (•7 -)n, 
the elements of which are the eigenvectors of the matrix DjM^D*. Let P' be the column 
matrix the elements of which are the components of the pressure pw in the new basis B' 
of Pjv(f^)- We have P' = R T P = R~ l P and the following equivalence 

(D^M~ l D)P = DjF A P' = R T DlF. (5.15) 

Since Z* N is a two dimensions space, B' 0 = (hi, h 2 ) is the basis of Z^ and (/i*) 3 <<<(JV+ 1)2 is 
a basis of Mjv = (Z^)^. Therefore, we determine the column vectors P' and P by 

p[ = p^ = o, Pi = ^—(R T D*F)i, 3<i<(N + l) 2 and P = RP', (5.16) 

and we have P^ s = P V.(r,s), 0 < r,s < N . In the same way as for the first discretization, 
we derive the column matrix U by the relation U = F — M~ l D*P, which gives, for 
0<j,k,t<N, 1 < u < N — 1, 




Ip 


- - £ PrW&Pr”, uff = Pin + ~£ PrU(r)P^- 

Pj r=0 P u r=0 
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2 N 2 N 

Finally, considering the boundary conditions, we set U t ’ 0 — U t ^ = 0, 0 < t < N, and we 
obtain the same formulas as in the first discretization for (divujv, divujv)iv, 

i 

(u — ujv, u — un)n and (p — pn,P — Pn)n • We also give the values of (divujv, divujv)jv, 

ii 

{p ~ PNiP ~ Pn)n an d ( u — u n, u — u n)n f° r N between 4 and 21. 


1 N 

4 

5 

6 

7 

8 

9 II 

T 

(divujv,divujv)jv 

1,84.10 -16 

5,45.HT 15 

3,60.10 16 

3,05.10 -15 

8,76.10- 16 

1,94.10 14 

1 

(. P~Pn,V-Pn)n 

0,246 

0,016 

0,019 

0,0039 

8 ,43.10~ 4 

1,55.10 3 

(u U/v, U U;y)jy 

0,62 

0,043 

0,054 

1,54. ur 3 

2,48.10~ 3 

3.54.10- 5 

N 

10 

11 

12 

13 

14 

15 

T 

(div ujv, divu w )^ 

1,24. ur 14 

9,28.10 15 

2,7.KT 13 

9,68.10“ 13 

4.25.10- 13 

3.43.10- 13 

1 

1 

2 ,3.10 5 

8,25.10~ 4 

4,33.10~ 7 

4,8.10~ 4 

5,9.ur 9 

3,l.HT 4 

1 

(U-Ujv,U-Ujv)a, 

T.0.10- 5 

5,7.10~ 7 

1,34.10 6 

6 ,79.10 -9 

1,85.10 8 

6,23.10~ n 

N 

16 

17 

18 

19 

20 

21 

T 

(divuAr.divuiv)^- 

4, 38.Hr 13 

2.24.10- 12 

3, 59.10 13 

1,9.10- 12 

8.71.10- 13 

(M 

7 

o 

oo 

CO 

of 

1 

(p-Pn,P-Pn)n 

6 ,11.10 11 

2,05.10~ 4 

4,96.10 13 

1,44.10 -4 

3,24.10~ 14 

1,04.10~ 4 

l 

(u-Ujv.U-Ujv)^ 

1,93.10~ 10 

4,53.HT 13 

1,57.10 12 

3.66.10- 14 

1 ,12.10- 13 

2,38.10 13 


Table 4.2 


i 

Logarithm to the basis 10 of (divujv, divujv)jy 
in function of N for the first and third discretizations 



Figure 4.3 
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CONCLUSION 

Continuous problem and discretizations 

We studied the Darcy problem by defining two equivalent variational formulations 
corresponding to two different bilinear forms b and b *. The first one requests a more regular 
pressure, which is bounded in H 1 ^). The second one is less classic and was introduced 
by A. Quarteroni and A. Valli (see [14]). This second variational formulation is important 
because it allowed us to constuct a spectral method which gives a divergence-free discrete 
solution. 

Moreover, by studying a mixed problem of Dirichlet-Neumann for the Laplace operator, 
we proved regularity results with the solution(u, p) in iL 1 (Q) 2 x H 2 (Q) as long as the data 
is regular enough. 

The first variational formulation led us to two spectral discretizations. In the first one, 
the discrete solution is not divergence-free, which is a disadvantage for some applications. 
However, we obtain a fully optimal error estimate for the velocity and the pressure. In 
the second discretization, it does not appear spurious modes for the pressure. Moreover, 
because of the exactness of the quadrature formulas, the error estimate is easier to obtain. 
However, the error estimate for the pressure is not optimal, because the inf-sup constant 
depends on N. 

The second variational formulation led to a third spectral discretization. There are 
spurious modes, which complicate the study, but the discrete velocity is divergence-free, 
which is important when the system is a stage of solving of a time-dependent problem, and 
the error estimate for the velocity and the pressure is fully optimal. In conclusion, this 
third discretization is the best discretization and we will only use it hereafter. 


Comparison of both spectral discretizations 

Tables 4.1 and 4.2 test the convergence of respectively the first and third discretiza¬ 
tions to the exact solutions. For the first discretization, we see the convergence to zero of 

i 

(div ujv, div ujvJat, about 10 -1 for N = 8 and about 10 -12 for N = 19. Concerning the 

i 

third discretization, we see that the quantity (div ujv, div uw)jv is very small for all values 
of N , about 10 -16 for N = 4 and about 10 -12 for N = 20. Thus, we verify that, in the 
third discretization, the discrete solution u^r is exactly divergence-free, which is important, 
as we saw previously. This property appears clearly in Figure 4.3, which represents the 

i 

logarithm to the basis 10 of the quantity (divujv,divu/v)jv as a function of A, for even N 
from 4 to 20, for both discretizations. 

Regarding the velocity, the discrete solution un converges fast to the exact solution u 
for both discretizations. For the pressure, we also obtain fast convergence, except for odd 
N in the third discretization. 
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